par(mfrow=c(2, 2))
data5_5 <- read.csv("data/习题数据（基于R，EXCEL格式）/csv/E5_5.csv")
ts_x <- ts(data5_5$pig_price)
# ts_x <- ts(data5_5$maize_price)

# 时序图
plot(ts_x, type="o", pch=5, col='#39CBB4')
# ADF检验，平稳性检验
# install.packages("aTSA")
library(aTSA)
adf.test(ts_x, nlag=3)

# 一阶差分
dif_x<-diff(ts_x)
plot(dif_x, type="o", pch=5, col='#39CBB4')
adf.test(dif_x, nlag=3)

# 白噪声，纯随机检验
for( k in 1:3) print(Box.test(dif_x, lag=6*k, type="Ljung-Box"))

# 自相关图
acf(dif_x, lag.max=20)
pacf(dif_x, lag.max=20)

# ARIMA((1, 2, 3, 14), 1, (1, 2, 3, 14, 15))
# 拟合梳系数模型 ARIMA((1, 2, 3, 14), 1, 0)
fit1 <- arima(ts_x, order=c(14, 1, 0), transform.pars = F, fixed = c(NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA))
fit1
fit1 <- arima(ts_x, order=c(14, 1, 0), transform.pars = F, fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA))
fit1

# ARIMA(0, 1, (1, 2, 3, 14, 15))
fit2 <- arima(ts_x, order=c(0, 1, 15), transform.pars = F, fixed = c(NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, NA))
fit2
fit2 <- arima(ts_x, order=c(0, 1, 15), transform.pars = F, fixed = c(0, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, NA))
fit2
fit2 <- arima(ts_x, order=c(0, 1, 14), transform.pars = F, fixed = c(0, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA))
fit2
fit2 <- arima(ts_x, order=c(0, 1, 14), transform.pars = F, fixed = c(0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA))
fit2

# 模型显著性检验
ts.diag(fit1)
ts.diag(fit2)
# 使用fit1预测
# 预测未来10期
library(forecast)
fore <- forecast(fit1, h = 10)
fore
plot(fore, lty=2)
lines(fore$fitted, col=4)
